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Abstract 

Hysteresis  in  smart  actuators  presents  a  challenge  in  control  of  these  actuators.  A  fundamental 
idea  to  cope  with  hysteresis  is  inverse  compensation.  But  due  to  the  open  loop  nature  of  inverse 
compensation,  its  performance  is  susceptible  to  model  uncertainties  and  to  errors  introduced  by 
inverse  schemes.  In  this  paper  we  develop  a  robust  control  framework  for  smart  actuators  by 
combining  inverse  control  with  the  l\  robust  control  theory.  We  show  that,  for  both  the  rate- 
independent  hysteresis  model  and  the  rate-dependent  one,  the  inversion  error  can  be  bounded  in 
magnitude  and  the  bound  is  quantifiable  in  terms  of  parameter  uncertainties  and  the  inversion 
scheme.  Hence  we  can  model  the  inversion  error  as  an  exogenous  disturbance  and  attenuate  its 
impact  by  robust  control  techniques.  Through  the  example  of  controlling  a  magnetostrictive  actua¬ 
tor,  we  present  a  systematic  controller  design  method  which  guarantees  robust  stability  and  robust 
trajectory  tracking  while  taking  actuator  saturation  into  account.  Simulation  and  experimental 
results  are  provided. 

Keywords :  Hysteresis;  Smart  Actuators;  Robust  control;  Inverse  compensation;  Actuator  satu¬ 
ration 

*  Corresponding  author. 
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1  Introduction 


Smart  materials,  such  as  magnetostrictives,  piezoelectrics,  shape  memory  alloys  (SMAs),  and  magne- 
torheological  (MR)  fluids,  all  display  coupling  phenomena  between  applied  electromagnetic/thermal 
fields  and  their  mechanical/rheological  properties.  Smart  actuators  and  sensors  made  of  these  mate¬ 
rials  have  been  receiving  tremendous  interest  due  to  their  broad  applications  in  areas  of  aerospace, 
manufacturing,  defense,  and  civil  infrastructure  systems,  to  name  a  few.  The  hysteretic  behavior 
widely  existing  in  smart  materials,  however,  makes  the  effective  use  of  these  actuators  and  sensors 
quite  challenging. 


Models  for  smart  actuators  that  capture  both  hysteresis  and  dynamic  behaviors  have  a  cascaded 
structure  as  shown  in  Figure  1(a)  [1],  where  W  is  a  hysteretic  operator  (with  possibly  some  other 
nonlinearities)  and  Ga( A)  is  a  linear  system.  In  this  paper  we  consider  the  discrete-time  setting  in  the 
interest  of  digital  control,  and  G( A)  denotes  the  A-transform  of  a  linear  time-invariant  (LTI)  system 
G.  We  recall  that  the  A-transform  (5(A)  is  just  the  usual  z-transform  of  G  with  A  =  z-1  [2]. 


fa) 


Controller 


(b) 


Figure  1:  (a)  The  model  structure  for  smart  actuators;  (b)  The  closed-loop  system  incorporating 
inverse  compensation. 

In  Figure  1(b),  Go  (A)  denotes  the  plant  to  be  controlled  by  the  actuator.  A  basic  approach  to  cope 
with  the  hysteresis  is  to  construct  an  (approximate)  right  inverse  operator  W~x  for  W,  then  u  «  u 
and  the  controller  design  problem  is  reduced  to  designing  a  linear  controller  K( A)  for  the  composite 
linear  system  Go  (A)  o  Ga(A).  The  idea  of  inverse  compensation  can  be  found  in,  e.g.,  [3,  4,  5,  6,  7]. 

The  most  popular  hysteresis  model  used  in  control  of  smart  actuators  has  been  the  Preisach 
operator  [3,  8,  9,  7].  The  Preisach  operator  provides  a  means  of  developing  phenomenological  models 
that  are  capable  of  producing  behaviors  similar  to  those  of  physical  systems.  For  a  detailed  treatment 
of  the  Preisach  operator,  we  refer  to  [10,  11,  12].  In  Part  I  of  this  two-part  paper  [13],  we  studied 
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modeling  and  inverse  control  of  smart  actuators  based  on  a  magnetostrictive  actuator.  In  particular, 
we  modeled  the  rateMndependent  hysteresis  with  essentially  a  Preisach  operator  alone,  and  the  rate- 
dependent  hysteresis  with  a  Preisach  operator  coupled  to  an  ordinary  differential  equation  (ODE). 
We  also  presented  identification  methods  and  efficient  real-time  inversion  schemes  for  these  models  in 

[13]. 

Due  to  the  open  loop  nature  of  inverse  compensation,  its  performance  is  susceptible  to  model 
uncertainties  and  to  errors  introduced  by  inversion  schemes.  To  combat  this  problem,  adaptive  in¬ 
verse  control  schemes  have  been  proposed  for  a  class  of  hysteresis  nonlinearities  with  parameterizable 
inverses  [4,  14,  15].  For  the  Preisach  operator-based  hysteresis  models,  however,  their  inverses  are  not 
parameterizable  in  general.  In  this  paper  we  develop  a  robust  control  framework  for  smart  actuators  by 
combining  inverse  control  with  the  l\  control  techniques.  We  show  that,  for  both  the  rate-independent 
hysteresis  model  and  the  rate-dependent  one,  the  inversion  error  can  be  bounded  in  magnitude  and 
the  bound  is  quantifiable  in  terms  of  parameter  uncertainties  and  the  inversion  scheme.  Hence  we 
can  model  the  inversion  error  as  an  exogenous  disturbance  and  attenuate  its  impact  by  robust  control 
techniques. 

The  design  requirements  for  the  controller  K( A)  can  be  roughly  stated  as:  in  the  presence  of  the 
inversion  error  and  the  uncertainties  in  Ga  and  Go,  for  all  desired  trajectories  in  a  certain  class, 

1.  The  closed-loop  system  is  stable; 

2.  The  tracking  error  is  minimized; 

3.  The  output  of  K  does  not  exceed  the  saturation  limits. 

We  take  the  saturation  constraint  (a  common  nonlinearity  in  actuators)  into  account  in  the  design  of 
K  to  ensure  that  the  overall  system  operates  in  the  linear  region  and  thus  predictions  based  on  the 
linear  design  are  credible.  The  controller  design  method  will  be  illustrated  through  the  example  of 
robust  trajectory  tracking  of  a  magnetostrictive  actuator. 

As  remarked  in  [13],  choosing  control  of  the  magnetostrictive  actuator  as  an  example  allows  us  to 
cover  both  the  rate-independent  hysteresis  case  and  the  rate-dependent  hysteresis  case.  The  approach 

presented  here  applies  directly  to  control  of  a  wide  class  of  smart  actuators. 

1  Throughout  the  paper,  the  word  “rate”  is  referred  to  how  fast  the  input  is  being  varied.  For  a  periodic  input,  the 
rate  is  directly  related  to  its  frequency. 
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The  remainder  of  the  paper  is  organized  as  follows.  In  Section  2  we  introduce  the  Preisach  operator 
and  review  an  identification  scheme  for  the  Preisach  operator.  In  Section  3  we  describe  the  model 
for  a  magnetostrictive  actuator.  We  discuss  quantification  of  bounds  on  inversion  errors  in  Section  4. 
We  then  formulate  and  solve  the  robust  control  problem  in  Section  5.  Simulation  and  experimental 
results  are  provided  in  Section  6.  Finally  we  conclude  in  Section  7. 

A  preliminary  version  of  some  results  in  this  paper  was  submitted  for  presentation  at  2003  American 
Control  Conference  [16]. 

2  The  Preisach  Operator 


For  a  pair  of  thresholds  (/?,  a)  with  (3  <  a,  consider  a  simple  hysteretic  element  7/3,cJ’>  •],  as  illustrated 
in  Figure  2.  For  u  G  C([0,T])  and  an  initial  configuration  £  G  {  —  1, 1},  the  function 

v  =  7/3,akC]  •  [0, T)  {-1,1} 


is  defined  as  follows  [11]: 


if  'u(O)  <  (3 
if  (3  <  u( 0)  <  a  , 
if  u( 0)  >  a 


and  for  t  G  (0,  T],  setting  Xt  =  {r  G  (0,  t\  :  u(r)  —  [3  or 


v(t)  ±  { 


«(0) 

-1 


if  Xt  =  0 

if  Xt  7^  0  and  'ii(maxlt)  =  (3 
if  Xt  ^  0  and  li(maxXt)  =  a 


v 


+1' 

_ _ 

1 

"t 

p 

1 

r 

a 

•  -1 

Figure  2:  The  elementary  Preisach  hysteron. 
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This  operator  is  sometimes  referred  to  as  an  elementary  Preisach  hysteron  (we  will  call  it  a  hysteron 
in  this  paper),  since  it  is  a  building  block  for  the  Preisach  operator. 

The  Preisach  operator  is  a  weighted  superposition  of  all  possible  hysterons.  Define 

V0  £  {(/?,«)  €R2  :(3<a}. 

Vo  is  called  the  Preisach  plane ,  and  each  (/?,  a)  G  Vo  is  identified  with  the  hysteron  73^.  For  u  G 
C([0,  T])  and  a  Borel  measurable  initial  configuration  (0  of  all  hysterons,  Co  :  Vo  — ►  {— 1, 1},  the  output 
of  the  Preisach  operator  Y  is  defined  as  [11]: 

z(t)  =T[u,Co](t)  =  [  'tf^a[u,Co(P,a)](t)du((3,a),  (1) 

Jv o 

where  v  is  a  finite,  signed  Borel  measure  on  Vo,  called  the  Preisach  measure. 

We  call  the  Preisach  measure  v  nonsingular  if  \u\  is  absolutely  continuous  with  respect  to  the 
two-dimensional  Lebesgue  measure.  By  the  Radon-Nikodym  theorem  [17],  if  v  is  nonsingular,  there 
exists  a  Borel  measurable  function  /i,  such  that 

r[«,Co \(t)=[  [  /J,(f3,a)%a[u,Co(f3,a)](t)dpda.  (2) 

J  JV o 

The  weighting  function  /a  is  often  referred  to  as  the  Preisach  function  [10]  or  the  density  function  [12]. 
We  consider  exclusively  the  case  fi>  0  in  this  paper. 

To  simplify  the  discussion,  throughout  the  paper  we  assume  that  /a  has  a  compact  support,  i.e., 
/i(/?,  a)  =  0  if  (3  <  (3o  or  a  >  ao  for  some  /?q,  ol o,  and  without  loss  of  generality,  we  let  ao  =  ~/3o  =:  ro  > 
0.  Then  it  suffices  to  consider  the  finite  triangular  area  V  =  {(/?,  <a)  G  M2|<a  >/?,/?>  —  vo,ol  <  ro}. 

At  time  £,  V  can  be  divided  into  two  regions:  V±(t)  =  {(/?,  a)  G  V\  output  of  7/3, a  at  t  is  ±1}.  In 
most  cases  of  interest,  each  of  V-  and  V+  is  a  connected  set  [10],  and  the  output  of  Y  is  determined  by 
the  boundary  between  V-  and  V+  if  the  Preisach  measure  is  nonsingular.  The  boundary  is  also  called 
the  memory  curve.  The  memory  curve  typically  has  a  staircase  structure  and  its  intersection  with  the 
line  a  =  f3  gives  the  current  input  value.  For  a  precise  characterization  of  the  set  T  of  memory  curves 
,  we  refer  to  [18].  The  memory  curve  ^0  at  t  =  0  is  called  the  initial  memory  curve  and  it  represents 
the  initial  condition  of  the  Preisach  operator. 

If  the  Preisach  measure  is  nonsingular,  we  can  identify  a  configuration  of  hysterons  with  a 
memory  curve  ^  in  the  following  way:  ^(/?,  a)  —  1  (— l,resp.)  if  (/3,a)  is  below  (above,  resp.)  the 


graph  of  rf>.  Note  that  it  does  not  matter  whether  takes  1  or  —1  on  the  graph  of  In  the  sequel 
we  will  put  the  initial  memory  curve  as  the  second  argument  of  T,  where  r[-,ipo\  =  r[-,^0]. 


A  constrained  least  squares  scheme  was  proposed  to  identify  the  Preisach  measure  v  in  [7].  In 
the  scheme,  the  input  is  discretized  into  L  +  1  levels  (called  discretization  of  level  L )  for  some  L  >  0 
and  that  leads  to  a  discretized  Preisach  operator  (Figure  3),  i.e.,  a  weighted  sum  of  finitely  many 
hysterons.  What  is  identified  in  [7],  is  a  collection  of  weighting  masses  sitting  at  centers  of  cells  in  the 
discretization  grid  (see  the  dark  dots  in  Figure  3).  We  can  then  obtain  a  nonsingular  approximation 
isp  to  the  true  Preisach  measure  v  by  assuming  that  each  identified  mass  is  distributed  uniformly  over 
the  corresponding  cell.  Note  that  the  density  fip  corresponding  to  vp  is  piecewise  uniform. 


Figure  3:  Discretization  of  the  Preisach  plane  (L  =  3). 


3  The  Model  for  Thin  Magnetostrictive  Actuators 


Magnetostriction  is  the  phenomenon  of  strong  coupling  between  magnetic  properties  and  mechanical 
properties  of  some  ferromagnetic  materials  (e.g.,  Terfenol-D).  Figure  4  shows  a  sectional  view  of  a  thin 
Terfenol-D  actuator.  By  varying  the  current  in  the  coil,  we  vary  the  magnetic  field  in  the  Terfenol-D 
rod  and  thus  control  the  displacement  of  the  rod  head.  Figure  5  displays  the  hysteresis  observed  in 
the  magnetostrictive  actuator. 

When  the  frequency  of  the  input  current  is  low  (typically  below  5  Hz),  the  magnetostrictive 
hysteresis  is  rate-independent  -  roughly  speaking,  the  shape  of  the  hysteresis  loop  is  independent  of 
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Figure  4:  Sectional  view  of  a  Terfenol-D  actuator  [19] (Original  source:  Etrema  Products,  Inc.). 


Figure  5:  Hysteresis  in  a  magnetostrictive  actuator, 
the  input  frequency,  and  it  can  be  modeled  by  essentially  a  Preisach  operator  2  [7]: 

mi)  =  com 

M(t)  =  r[H(-),iM(t)  ,  (3) 

y(t)  =  cMM2(t) 

where  /  is  the  input  current,  y  is  the  displacement  of  the  actuator  head,  M  and  H  are  the  bulk 
magnetization  and  the  magnetic  field  (assumed  uniform)  along  the  rod  direction,  respectively,  T  is  the 
Preisach  operator,  and  Co  and  cm  are  positive  constants. 

When  the  input  frequency  gets  high,  the  magnetostrictive  hysteresis  is  rate-dependent.  Venkatara- 
man  and  Krishnaprasad  proposed  a  bulk  magnetostrictive  hysteresis  model  for  the  thin  rod  actuator 
based  on  energy  balancing  principles  [20,  19].  The  model  has  a  cascaded  structure  as  shown  in  Fig¬ 
ure  6.  Note  the  resemblance  of  Figure  6  with  Figure  1(a).  W  takes  care  of  the  M  -  H  hysteresis  and 

2To  ease  discussion,  we  have  absorbed  the  (constant)  bias  field  Hbias  [13]  into  I(t). 


7 


the  eddy  current  losses,  and  the  magnetoelastic  dynamics  of  the  rod  is  lumped  into  a  second  order 
linear  system  G(s).  G(s)  has  a  state  space  representation  [20,  19] (after  some  manipulations): 


y(t)  +  2&0y(t)  +  ioly(t)  =  iolcMM2(t ), 
where  uj o  and  £  are  positive  constants. 


(4) 


Figure  6:  Model  structure  of  a  magnetostrictive  actuator. 

By  replacing  the  switching  ODE  model  in  [20,  19]  with  a  Preisach  operator  T  for  the  M  -  H 
hysteresis,  we  have  proposed  a  new  dynamic  model  [13,  18]  for  the  W  block: 

f  H(t)  +  M(t)  =  c1(m-Z&) 

where  c\  is  a  positive  constant.  We  note  that  the  Preisach  operator  is  coupled  to  an  ODE  in  an 
unusual  way  in  (5). 

We  also  note  that  setting  derivatives  in  (4)  and  (5)  to  zero,  the  dynamic  model  degenerates  to  the 
rate-independent  hysteresis  model  (3).  More  details  on  magnetostrictive  actuators  and  the  modeling 
can  be  found  in  [13]. 

4  Quantification  of  the  Inversion  Error 

Recall  Figure  1(b).  In  general  W~l  is  not  an  exact  (right)  inverse  of  W,  and  two  factors  may  contribute 
to  the  inversion  error  eu:  parameter  uncertainties  in  W  and  non-existence  of  exact  inverse  schemes. 

There  are  two  possible  ways  to  model  eu.  One  is  to  model  it  as  the  output  of  some  uncertainty 
block  A  (Figure  7(a)),  and  the  other  is  to  simply  model  it  as  an  exogenous  disturbance  v  (Figure  7(b)). 
For  the  Preisach  operator-based  models,  eu  is  independent  of  u  and  it  is  possible  that  eu  ^  0  for  u  =  0. 
Therefore  there  exists  no  stable  A  such  that  eu  =  A u,  and  we  will  treat  eu  =  v  as  an  external  noise. 

We  need  specify  the  signal  spaces  for  quantification  of  the  inversion  error.  The  inversion  error  for 
the  Preisach  operator  is  bounded  in  magnitude  instead  of  in  energy.  Hence  a  natural  choice  for  the 


Figure  8:  The  error  in  inversion  of  the  Preisach  operator. 

signal  spaces  is  and  not  Z2.  Also  it  is  more  appropriate  to  use  Z^  for  the  desired  trajectory  and  the 
tracking  error.  Another  advantage  of  using  l ^  for  signals  is  that  the  actuator  saturation  constraint  can 
be  easily  handled  in  the  corresponding  l\  robust  control  theory,  while  it’s  very  hard  to  be  formulated 
in  Hqo  control  theory. 

We  now  quantify  the  error  bounds  in  inversion  of  the  Preisach  operator  and  the  dynamic  model 
(5).  Here  we  are  concerned  with  cm[-\  —  M[-]  —  M[*],  where  M[-\  and  M[-\  denote  the  trajectories  of 
achieved  output  and  desired  output  of  the  Preisach  operator,  respectively.  The  bound  on  eu  when  the 
square  nonlinearity  in  Figure  6  is  included  can  be  easily  derived  from  the  bound  on  e^. 

We  will  not  review  the  inversion  algorithms  here,  since  they  have  been  given  in  [13].  Interested 
readers  are  referred  to  [13]  or  [18]  for  details. 


4.1  Error  in  inversion  of  the  Preisach  operator  with  a  nonsingular  measure 

Consider  Figure  8,  where  V  is  a  Preisach  operator  with  nonsingular  Preisach  measure  v. 

Given  a  sequence  {M[fc]}  G  ZQ c  and  the  initial  memory  curve  if;  oG$  (the  set  of  memory  curves), 
the  inversion  problem  for  the  Preisach  operator  is  to  find  the  input  sequence  {H[k]}  G  ZQ c,  such  that 
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the  output  sequence  {M[k]}  matches  {M[/c]},  i.e., 

M[k]=T[H[]^0][k}. 

Remark  4.1  Note  that  we  now  have  a  sequence  instead  of  a  continuous  time  function  as  the  first 
argument  of  T.  To  avoid  confusion ,  we  tacitly  understand  that  the  input  is  monotonically  changed 
from  H[k\  to  H[k  +  1],  Throughout  the  paper  we  may  use  a  sequence  or  a  continuous  time  function 
as  the  first  argument  of  T  depending  on  the  context. 

As  discussed  in  [13],  for  a  Preisach  operator  with  a  nonnegative,  nonsingular  Preisach  measure,  the 
inversion  problem  can  be  solved  with  arbitrary  accuracy  if  the  Preisach  measure  is  known  exactly.  In 
this  case,  it  suffices  to  consider  the  inversion  problem  of  length  1:  given  G  and  M  G  [— Ms,  Ms], 
where  Ms  is  the  saturation  output  of  the  Preisach  operator,  find  H  such  that 

M  =  T[H,  V’o].  (6) 


4.1.1  Error  due  to  the  inversion  scheme 

If  the  Preisach  measure  v  is  known  exactly  but  otherwise  has  a  general  distribution,  an  iterative  algo¬ 
rithm  is  available  to  solve  (6)  with  arbitrary  accuracy  [13].  Let  the  stopping  criterion  be  <  e, 

where  M ^  denotes  the  output  value  of  T  achieved  at  the  n-th  iteration.  Then  it’s  straightforward 
that 

||  eM  ||oo<  e  for  any  M[\  G  l^. 


4.1.2  Error  due  to  the  parameter  uncertainty 

If  the  Preisach  measure  v  is  unknown,  we  can  obtain  a  nonsingular  approximation  vp  with  a  piecewise 
uniform  density  fip  as  discussed  in  Section  2.  The  Preisach  operator  with  measure  vp  can  be  inverted 
exactly  (in  a  finite  number  of  steps)  [13].  Hence  the  inversion  error  eM  is  now  solely  due  to  the 
measure  uncertainty  \v  —  vv\.  It  turns  out  that  we  can  quantify  the  error  bound  in  terms  of  the 
relative  identification  error  and  the  discretization  level  L  of  the  Preisach  plane: 

Proposition  4.1  Let  the  true  Preisach  measure  v  be  nonnegative  and  nonsingular  with  density  pi. 
Let  fi  be  bounded  by  some  constant  jd  >  0.  Given  a  discretization  of  level  L,  denote  the  integral  of  /a 
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over  a  cell  i  as  vf,  1  <  i  <  Nc,  where  Nc  is  the  total  number  of  cells.  Denote  by  ^  the  identified  mass 
for  cell  i.  Assume  the  relative  error  in  identification  is  5j,  i.e .,  z  0  *  <5/,  1  <  i  <  Nc.  Then 

^ i 

||  CM  ||oo<  SiMs  + 

where  Ms  is  the  saturation  output  of  the  Preisach  operator  with  measure  v ,  and  ro  is  as  defined  in 
Section  2. 

Proof.  Define  fiv  as  discussed  earlier.  We  obtain  another  Preisach  measure  with  a  piecewise  constant 
density  /i®  by  distributing  vf  uniformly  over  the  cell  i,  1  <  i  <  Nc.  To  distinguish  the  Preisach 
operators,  we  will  put  the  corresponding  density  as  the  subscript  of  T,  e.g.,  T ^  denotes  the  Preisach 
operator  with  Preisach  density  /a. 

Given  M[-\  G  l oo  and  we  denote  the  output  of  T_1  in  Figure  8  as  H,  where  T_1  is  the  inverse 
of  rMp.  Then,  Vk  >  0, 

l^[fc]|  =  |rM[^^o][fc]  -r^[^^o][fc]| 

<  \r,[H^0}[k\-r^[H^0}[k\\  +  \T^[H^][k\  -T^[H,M[k\\.  (7) 

All  three  Preisach  operators  involved  in  (7)  share  the  same  memory  curve  fi[k\,  VA:  >  0.  It’s  obvious 

that  the  second  term  of  (7)  is  bounded  by  5jMs.  To  bound  the  first  term,  we  note  that  for  any  k  >  0, 

the  memory  curve  fi[k\  spans  L  —  1  square  cells  and  one  triangular  cell  (Figure  9).  Any  cell  not  touched 

by  fi[k\  will  contribute  the  same  amount  to  T ^[H,  fi>o\[k\  and  T^o  [77,  V’o]^]-  Hence  the  first  term  of  (7) 

is  bounded  by  twice  the  interal  of  /a  over  cells  spanned  by  fi[k\,  which  is  further  bounded  by 

2/2(2r0)2(L  ~\)  ^  8 p,r% 

L2  L  ' 

This  completes  the  proof.  □ 

Remark  4.2  From  Proposition  fil,  the  bound  on  the  inversion  error  consists  of  two  parts:  the  first 
part  is  proportional  to  the  relative  identification  error ,  and  the  second  part  is  inversely  proportional  to 
the  level  L  of  discretization. 

Remark  4.3  The  assumption  that  fio  is  known  is  not  very  restrictive  since  in  many  cases  we  have 
the  choice  to  initialize  the  system.  On  the  other  hand,  if  o  is  not  known  exactly,  we  can  easily  include 
a  term  in  ||  ||oo  which  takes  care  of  the  uncertainty  in  fio. 
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Figure  9:  Illustration  of  the  proof  of  Proposition  4.1  (L  —  8). 

4.2  Error  in  inversion  of  the  dynamic  hysteresis  model 


Given  the  initial  memory  curve  i/jq  and  a  trajectory  of  M(-),  inversion  of  (5)  is  to  find  /(•),  so  that 
the  output  of  the  block  W  in  Figure  6  is  M(-).  The  following  (formal)  inversion  scheme  was  proposed 
in  [13]: 

I{t)  =  ±-(H(t)  +  M(t))  +  ^,  (8) 

Cl  Co 

where  H(t)  —  r_1[M(-),  ^o](t).  But  if  there  is  uncertainty  in  the  model  parameters,  it  is  very  hard  to 
derive  a  bound  for  the  inversion  error.  We  now  present  another  inversion  algorithm.  This  algorithm 
leads  to  an  inversion  error  even  if  the  exact  parameters  are  known,  but  it  will  allow  us  to  quantify  the 
inversion  error  when  model  uncertainties  are  present. 


Eq.  (5)  can  be  rewritten  as: 

{  m  =  _  (9) 

where  g(t)  carries  the  interpretation  of  at  time  t,  and  it  depends  on  both  the  state  ^  (the  memory 

curve  at  t)  and  the  sign  of  H  [18].  When  the  Preisach  measure  is  nonnegative  and  nonsingular  with  a 
piecewise  continuous  density,  we  have  0  <  g(t)  <  Cgi  for  some  constant  Cg  >  0  [18].  We  can  view  (9) 
as  perturbed  from  the  following  decoupled  system: 


m  = 

M(i)  =  r[H(.),V>o](t) 


where  g  G  [0 ,Cg\  is  some  constant  (at  our  disposal).  Based  on  (10),  an  approximate  inversion  scheme 
for  (9)  is  given  formally  by 


m  =  ^ H(t )  +  ^ 


(ii) 
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We  have  two  ways  to  implement  the  second  equation  in  (11)  in  the  discrete  time,  which  correspond 
to  the  explicit  Euler  scheme  and  the  implicit  Euler  scheme  in  discretizing  the  first  equation  in  (10), 
respectively:  for  k  >  0, 

m  =  H[k }  -  H[k  -  1])  +  (12) 

c\h  co 

m  =  -  s[k  - 1])  +  (is) 

ci  ti  Co 

where  h  is  the  time  step  size,  H[—  1]  =  H[ 0]. 

Remark  4.4  Direct  discretization  of  the  first  equation  in  (10)  by  the  explicit  Euler  scheme  (a  similar 
remark  applies  to  the  implicit  Euler  scheme)  gives: 

m  =  +  1]  -  H[k])  +  EH , 

C\tl  Co 

but  this  is  not  a  causal  system  and  thus  not  realizable.  An  intrinsic  delay  is  introduced  in  the  inversion 
due  to  the  dynamics  in  the  rate- dependent  hysteresis  model. 

We  now  want  to  study  the  errors  caused  exclusively  by  the  inversion  algorithms,  i.e.,  we  assume 
that  we  have  exact  values  of  parameters.  For  the  algorithm  (12),  the  discrete-time  version  of  the  first 
equation  in  (9)  is  obtained  by  the  explicit  Euler  scheme: 

=  «  (,w- ffi,.  (i4) 

h  1  +  g[k\  co 

where  g[k\  =  g{kh).  Similarly,  if  the  inversion  algorithm  (13)  is  used,  we  will  use  the  corresponding 
discrete-time  model  obtained  by  the  implicit  Euler  scheme.  Figure  10(a)  shows  the  problem  setup  for 
the  explicit  Euler  case. 

For  the  purpose  of  deriving  the  bound  on  the  inversion  error,  we  will  not  need  the  exact  values  of 
g[k\.  Due  to  the  delay  caused  by  the  inversion,  the  error  is  now  defined  as  (Figure  10(d)): 

eM[k]  =  M[k\  -  M[k  -  1]. 


Proposition  4.2  Let  the  Preisach  measure  be  nonnegative  and  nonsingular  with  a  piecewise  con¬ 
tinuous  density  (i.  Let  the  Preisach  operator  V  be  Lipschitz  continuous  with  Lipschitz  constant  L^. 
Consider  the  inversion  algorithm  obtained  from  the  explicit  Euler  method  (Figure  10(a)).  Let  H[—  1]  = 
H[ 0]  =  ^[0].  Pick  g  G  [0,  Cg\.  Then  for  any  M[f]  G  Zoo,  for  any  fio  £ 

||  1 1 oo  <  2LM7er0,  (15) 
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(b) 


(c) 


M[k] 


M[k-1]  X  M[k]  ^ 


(d) 


Figure  10:  The  error  in  inversion  of  the  rate-dependent  hysteresis  model. 


where 


ma x{g,  rM  } 


7e  = 


-  -  c 

The  optimal  g  to  minimize  ye  is  c 


1  —  maxj^1  —  1,1  — 

L  CO 


hci 

co(l+Cg) 


} 


Proof.  We  first  derive  a  bound  for  e#,  defined  by  e#[fc]  =  H[k]  —  H[k  —  1],  k  >  0.  Substituting  (12) 
into  (14),  we  have 

eif[fe  +  1]  =  a[fc]e#[fc]  +  b[k\(H[k\  -  H[k  -  1]),  (16) 


where 


From  (16),  we  compute 


a[k]  =  1  - 


hc\ 


co(i  +  g[k\Y 


b[k]  ±  9 


9[k] 


1+#]' 


eH[k  +  1] 


k 


k  k 


(II«[i])«[0]  +  E(  II  «b]»l(tf[i]  -  H[i  -  1]). 

2=0  2=0 


(17) 
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Since  e#[0]  =  0, 


\en[k  +  1]  | 


<  2('^ai)b  ||  H  Hoc 

i= 0 

2b 


< 


1  —  a 


H  ||oo? 


(18) 


where 


It’s  easy  to  verify  that 


a  =  max  1 1 

xe[o  ,cg\ 


hc\ 

C0(l  +  X ) 


b  —  max 

xe[o  ,cg\ 


g-x 

\  X 


Ml  1  hcx  -  Cg-g 

a  =  max{— - 1, 1  -  ~j— b  =  max{g,  }. 

Co  Co(l  +  Cg)  1  +  Cg 

Therefore  ||  en  ||oo^  2ye  ||  H  ||oo.  The  error  en  can  be  thought  of  as  the  output  of  some  uncertainty 
block  A h  with  the  induced  gain  less  than  or  equal  to  2ye  (Figure  10(b)).  But  since  T,T_1  sit  outside 
the  dashed  box  in  Figure  10(b),  we  can  not  carry  Ah  along  further.  Instead  we  represent  e#  as  an 
exogenous  disturbance  with  magnitude  bounded  by  27ero  (Figure  10(c)).  Eq.  (15)  now  follows  using 

the  Lipschitz  continuity  and  the  time  invariance  properties  of  T.  It’s  easy  to  see  that  the  optimal  g 

c 

minimizing  the  error  bound  is  c  9+2.  □ 


Similarly  we  can  derive  the  error  bound  for  the  implicit  Euler  algorithm  (13): 


Proposition  4.3  Let  the  assumptions  in  Proposition  2^.2  hold.  Consider  the  implicit  Euler  algorithm 
(13).  Then  for  any  M[f]  G  lQ 0,  for  any 

||  | loo  <  2LM7jr0,  (19) 


where 


7 i  =  max{ 


Cg  —  g  .  cq(1  +  Cg)  +  c\h 


1  + 


c\h 

co 


1  +  CQ  + 


c\h 

co 


c\h 


The  optimal  g  to  minimize  7 i  is 


(CQ+Clh)Cg 

2(co+cih)+coCg  * 


Remark  4.5  For  the  explicit  algorithm,  the  step  size  h  has  to  be  chosen  small  enough  to  ensure 
stability  of  (If)  and  (16).  The  implicit  algorithm,  however,  is  stable  Mh  >  0.  Therefore  the  implicit 
algorithm  is  preferred  in  general. 


Remark  4.6  Propositions  f.2  and  f.3  use  the  Lipschitz  continuity  of  the  Preisach  operator.  We 
recall  from  [11]  that  if  the  Preisach  measure  is  nonsingular,  the  Preisach  operator  will  be  continuous. 
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Furthermore,  if  the  measure  v  satisfies  the  condition 


v(N(iP,e))<eL^\/^e*,  (20) 

for  some  L y  >  0,  where  N(fi>,e)  denotes  the  e-neighbourhood  of  ip,  then  the  Preisach  operator  is  Lip- 
schitz  continuous  with  the  Lipschitz  constant  2 Ly.  The  condition  (20)  is  usually  satisfied  in  practice. 


Propositions  4.2  and  4.3  quantify  the  errors  solely  due  to  inversion  algorithms.  It’s  straightforward  to 
extend  the  error  estimates  to  the  case  that  there  are  parametric  uncertainties  in  Co  and  c\.  The  error 
due  to  inversion  of  the  Preisach  operator  and  to  the  uncertainty  in  the  Preisach  measure  can  also  be 
included  as  done  in  Subsection  4.1. 

When  the  square  operator  is  present,  like  in  the  case  of  a  magnetostrictive  actuator,  the  estimate 
of  eu  can  be  derived  from  that  of  e^.  Let  u  E  (recall  Figure  1(b)),  with  umin  >  0.  One 

can  easily  verify  that 

||  &u  Hoo^ll  &M  ||oo  T  ^  II  ||oo  \J  ^ max  • 

5  Formulation  of  the  Robust  Control  Problem 


In  this  paper,  we  consider  Go  (A)  to  be  the  identity  operator,  i.e.,  we  are  interested  in  trajectory  tracking 
of  the  actuator  head  itself.  Figure  11  shows  the  closed- loop  system  after  the  inverse  compensation  is 
done,  where  the  exogenous  noise  v  represents  the  inversion  error.  From  the  previous  section,  ||  v  ||oo<  v, 
and  v  is  quantifiable  in  terms  of  inverse  schemes  and  parametric  uncertainties.  Ga(A)  stands  for  the 
discretized  version  of  G(s )  in  Figure  6.  The  composition  A  o  Wo  (A)  represents  the  deviation  of  the 
actual  plant  from  the  nominal  plant  Ga(A).  We  assume  that  A  can  be  any  nonlinear  operator  with 
||  A  || iOQ-ind<  where  ||  •  || i^-ind  denotes  the  induced  operator  norm  when  the  signal  space  is  l0 c. 
Wo  (A)  is  a  weighting  function  and  it  reflects  that  the  model  uncertainty  is  larger  at  a  higher  frequency. 

Let  ||  yref  ||oo<  r,  where  yref  is  the  reference  trajectory.  The  error  ey  =  yref  —  y  is  fed  into  the 
controller  K{ A).  The  delay  A  following  K{ A)  is  due  to  inversion  of  the  dynamic  hysteresis  model. 
Let  the  saturation  limits  of  the  actuator  be  —u  and  u  respectively.  Then  the  saturation  constraint 
translates  into  ||  uq  ||oo^  1,  where  uo  is  as  defined  in  Figure  11.  The  case  umin  —umax  will  be 
discussed  in  Section  6. 
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Figure  11:  Robust  control  of  a  magnetostrictive  actuator. 


Figure  12:  Formulation  of  the  robust  control  problem. 


There  are  two  delays  in  the  loop  since  Ga( A)  contains  a  pure  delay.  This  motivates  us  to  define 
the  tracking  error  eo  as 

e„[k}  = 

7 

where  7  >  0  is  the  desired  disturbance  attenuation  level.  To  ease  the  formulation,  we  normalize  signals 


v  and  yrefi  and  regard  vq  and  ro  as  inputs  to  the  system  with  ||  ^0  ||oo<  1?  ||  ^0  ||oo<  1  (Figure  11). 


GW 


In  terms  of  G,  the  closed-loop  system  in  Figure  11  can  be  simplified  as  in  Figure  12  (a). 


The  control  objective  is:  find  the  smallest  7  and  a  stabilizing  controller  K( A),  such  that 
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1.  the  closed-loop  system  is  stable  for  any  A  with  ||  A  \\iOQ-ind<  1, 

2.  ||  eo  ||oo^  1  if  A  =  0,  V^o,  ro  with  ||  vq  ||oo<  1  and  ||  ro  ||oo<  1,  and 

3.  ||  u0  ||oo<  1  if  A  =  0,  Vr0,  r0  with  ||  r0  ||oo<  1  and  ||  r0  ||oo<  1- 

If  we  define  the  exogenous  input  w  and  the  regulated  output  z  as 


then  items  2  and  3  above  are  equivalent  to  ||  &zw  ||i<  1,,  where  <$>zw  denotes  the  transfer  function 
from  w  to  z,  and  ||  •  ||i  denotes  the  li  norm  of  a  LTI  system  [2].  By  the  small  gain  theorem,  (5)  is 
equivalent  to  requiring  robust  stability  of  the  system  when  we  wrap  a  nonlinear  uncertainty  block  A p 
from  z  to  w  with  ||  A p  Wi^-ind  <  1,  as  shown  in  Figure  12  (b). 

Now  the  control  problem  can  be  reformulated  as:  find  the  smallest  7  and  a  stabilizing  controller 
K( A),  such  that  the  closed- loop  system  in  Figure  12  (b)  is  robustly  stable  for  all  A  G  A,  where 
A  =  {A  =  diag(A,  Ap)  :  A  is  nonlinear  and  of  dimension  1  X  1,  Ap  is  nonlinear  and  of  dimension 
2x2,  ||  A  H/^— 1}. 

To  solve  the  robust  control  problem,  we  need  determine,  for  a  fixed  7  >  0,  whether  we  can  find 
a  stabilizing  K( A),  such  that  the  closed-loop  system  is  stable  for  all  A  G  A.  This  will  be  called  the 
robust  control  problem  with  disturbance  attenuation  level  7,  and  it  is  solvable  if  and  only  if 

inf  _  inf  ||  D~1Fi(G,K)D  ||i<  1,  (23) 

stabilizing  k  D 

where  D  =  {D  —  diag(d\,  o?2?  ^2)  •  ^1,^2  >  0},  and  denotes  the  lower  Linear  Fractional 

Transformation  [2].  We  now  give  a  sketch  of  how  to  solve  (23).  Interested  readers  are  referred  to  [2] 
for  details  of  the  approach. 

We  restrict  ourselves  to  finite  dimensional  LTI  (FDLTI)  controllers.  Simultaneous  optimization 
with  respect  to  D  and  K  is  hard,  and  a  practical  approach  is  to  use  the  D  —  K  iteration  method  and 
decompose  the  problem  (23)  into  a  sequence  of  decoupled  optimization  problems. 

In  Step  1  of  each  D  —  K  iteration,  for  a  fixed  Dg  D,  we  want  to  solve 

inf  ||  D~1Ft(G,K)D  ||i  .  (24) 

stabilizing  k 
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Partition  G  into  a  2  x  2  block  matrix  as  shown  in  (22)  and  denote  it  as 

6Jgu  g12\ 

\  G‘21  G22  ) 

Since  G22  for  our  problem  is  stable,  the  set  of  stabilizing  FDLTI  controllers  K  is  parametrized  by: 

K( A)  =  S— ,  Q  G  nnlZ\  (25) 

1  —  G22  Q 

and  the  scaled  achievable  closed-loop  maps  is  parametrized  by 

D-1F/(G,  K)D  =  E  -  UQV,  Q  €  nn^\  (26) 

where  E  =  D~1GnD,  U  =  D~lG\2,  V  =  G21D ,  and  TZH^1  denotes  the  space  of  real  rational 
(complex)  functions  that  are  analytic  on  the  open  unit  disc  and  bounded  on  the  unit  circle.  Therefore 
(23)  is  transformed  into 

inf  \\E-UQV\h.  (27) 

Qenn1^1 

Problem  (27)  is  a  multi-block  li  model  matching  problem  and  we  can  approximate  it  by  a  one-block  li 
model  matching  problem  through  delay  augmentation  (DA).  The  latter  problem  is  then  solved  using 
linear  programming.  We  also  obtain  a  sub-optimal  controller  for  (26)  from  the  DA  method. 


Remark  5.1  For  the  system  we  consider,  the  only  zeros  that  the  delay  augmented  matrices  Ujy  and 
VN  (c.f.  [2])  have  inside  the  unit  disk  are  0’s.  This  has  two  pleasant  consequences: 

1.  In  computation  of  null  chains  and  evaluation  of  the  zero  interpolation  conditions,  relevant  coef¬ 
ficients  can  be  obtained  directly  from  the  impulse  responses  and  we  thus  avoid  expensive  symbolic 
calculation  of  high  order  derivatives. 

2.  From  the  zero  interpolation  conditions,  the  upper  bound  on  the  (finite)  length  of  the  closed-loop 
impulse  response  can  be  derived  explicitly. 


In  Step  2  of  each  D  —  K  iteration,  an  analytical  expresssion  for  the  optimal  D*  exists  since  there 
are  only  two  blocks  in  the  structured  uncertainty  class  A. 
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Figure  13:  Effect  of  the  model  uncertainty  on  7*. 

6  Simulation  and  Experimental  Results 


6.1  Effects  of  design  parameters  on  the  optimal  disturbance  attenuation  level 


We  first  present  some  computation  results  on  how  the  optimal  attenuation  level  7*  is  affected  by 
the  following  factors:  the  magnitude  of  uncertainty,  the  magnitude  v  of  the  inversion  error,  and  the 
saturation  limit  u. 


The  sampling  frequency  we  use  is  2000  Hz.  The  corresponding  Ga( A)  is 

*  _  2.23  x  10-11  A2  +  4.28  x  10-nA 

Ga(A)  ”  0.147A2  —  0.549A  +  1  ' 

We  choose  the  continuous  time  weighting  function  to  be  Wq(s)  =  7^^,  where  cw  >  0  determines 
the  magnitude  of  the  uncertainty  in  the  plant.  Discretizing  Wq(s)  gives 


Wb(A) 


1.1759c™  (A  -  1.0005) 
A  -  1.1765 


We  let  r  30. 


Figure  13  shows  the  effect  of  the  uncertainty  magnitude  on  7*.  Other  parameters  used  are  v  = 
0.1MS2,  u  —  7.5M2,  where  Ms  is  the  saturation  magnetization.  Since  the  range  of  u  in  the  case  of 
a  magnetostrictive  actuator  is  [0,  M2]  (recall  Figure  6),  expressing  v  and  u  in  terms  of  M2  allows 
one  to  make  more  concrete  sense  out  of  these  numbers.  From  Figure  13,  we  see  that  the  higher  the 
uncertainty,  the  bigger  7*. 
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Figure  14:  Effect  of  the  inversion  error  on  7*. 


Figure  15:  Effect  of  the  saturation  limit  on  7*. 

Figure  14  displays  how  7*  varies  with  the  magnitude  v  of  the  inversion  error,  where  we  have  fixed 
cw  =  6.53  x  10-13  and  u  =  1.25 Mg.  As  one  expects,  the  optimal  attenuation  level  7*  increases  as  v 
increases. 

Figure  15  shows  how  7*  is  affected  by  the  saturation  constraint.  We  have  used  cw  =  6.53  x  10-13 
and  v  0.1  Mg.  7*  drops  when  u  increases,  but  7*  becomes  a  constant  when  u  hits  4.5M^,  beyond 
which  the  saturation  constraint  no  longer  plays  a  role. 
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6.2  Results  of  trajectory  tracking 


The  saturation  constraint  considered  so  far  is  of  the  form  \u\  <  u.  But  for  real  actuators,  the  saturation 
limits  may  be  asymmetric,  i.e.,  umin  ^  For  example,  for  magnetostrictive  actuators,  u  G 

[0,  M2].  To  handle  the  general  constraint  u  G  [umin,  umax\,  we  let  u  =  Umax~Umin  and  u 5  =  Umax+Umin  . 
The  quantity  u  is  the  saturation  limit  to  be  used  in  the  controller  design,  while  u 5  is  a  bias  input  to 
be  injected  into  the  system.  Then  the  actual  control  will  be  u  =  uc  +  with  \uc\  <  u. 


Since  the  gain  of  Wq  is  close  to  0  for  a  dc  signal,  we  can  ignore  the  contribution  of  u 5  to  the  actuator 
output  y  through  the  A  o  W0  branch.  Its  contribution  through  the  Ga  branch  can  be  calculated  as 

AGa(A) 


Vb  = 


-ub. 


1  +  XK(X)Ga(X) 

The  previous  robust  control  framework  applies  if  we  add  yi  to  the  reference  trajectory  yref  (or  alter¬ 
natively,  taking  off  from  y). 


As  we  have  seen  from  Figure  15,  the  tracking  performance  deteriorates  as  the  saturation  constraint 
u  is  tightened.  For  the  magnetostrictive  actuator,  u  =  0.5 M2  and  strictly  enforcing  this  constraint 
will  lead  to  large  tracking  errors.  This  reveals  the  limitation  of  pure  linear  design  for  an  intrinsically 
nonlinear  plant.  Hence  a  practical  approach  would  be  to  properly  relax  the  constraint. 


Figures  16(a)  and  17(a)  show  the  simulation  results  of  tracking  a  sinusoidal  signal  and  an  irregular 
signal,  respectively.  The  current  /  applied  is  also  displayed.  The  controller  K( A)  is  designed  based 
on  f  =  30,  cw  =  3.3  x  10-13,  v  =  0.1  M^,  and  u  =  3.25 M%.  Figure  18  shows  the  output  of  K( A)  in 
tracking  the  sinusoidal  signal,  and  we  see  that  although  we  set  u  —  3.25 in  the  controller  design, 
the  control  stays  in  the  (true)  unsaturated  region  [— 0.5M^,  0.5M^]  except  during  the  transient  period 
at  the  beginning. 


Our  composite  controller  (the  linear  robust  controller  plus  the  inverse  algorithm)  is  compuation 
efficient  and  we  can  implement  it  in  real-time.  Figures  16(b)  and  17(b)  show  the  experimental  results 
of  trajectory  tracking  based  on  the  same  controller  used  in  the  simulation.  It  matches  well  with  the 
simulation  result  and  the  overall  performance  is  satisfactory. 


The  saturation  limit  u  can  not  be  “over-relaxed” .  For  example,  we  design  another  controller  based 
on  f  —  25,  cw  =  3.3  x  10— 13 ,  v  =  0.05M^,  and  u  =  5MS2.  The  simulation  result  (Figures  19(a)) 
based  on  this  new  controller  is  better  than  that  in  Figure  16(a).  But  when  we  put  the  controller  into 
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(a) 


(b) 


Figure  16:  Results  of  tracking  a  sinusoidal  signal. 


(a)  Simulation  result;  (b)  Experimental  result. 


(a) 


(b) 


Figure  17:  Results  of  tracking  an  irregular  signal. 


(a)  Simulation  result;  (b)  Experimental  result. 
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Figure  18:  The  control  output  uq  in  simulation  of  tracking  the  sinusoidal  signal. 

the  experiment,  the  tracking  performance  suffers  from  the  persistant  saturation  (Figure  19(b)).  This 
justifies  the  necessity  of  including  the  saturation  constraint  in  the  problem  formulation. 


- Desired  trajectory 

J -  Measured  trajectory 

r  \  :  \  \  j 

. \Jj :  \Ji  . \Jf 

. . y  vf . y  y 

o1 - 1 - 1 - 1 - 1 - 1 - 1 - 1 - 1 - 1 - 

0  0.05  0.1  0.15  0.2  0.25  0.3  0.35  0.4  0.45  0.5 

Time  (sec.) 


Time  (sec.) 


(a)  (b) 

Figure  19:  Results  based  on  an  “over-relaxed”  controller,  (a)  Simulation  result  of  trajectory  tracking; 
(b)  Experimental  result  of  trajectory  tracking. 


7  Conclusions 


In  this  paper,  we  presented  a  robust  control  framework  for  smart  actuators  by  combining  the  inverse 
compensation  with  the  linear  robust  control  theory,  where  the  essential  ideas  were  illustrated  through 
the  example  of  robust  tracking  of  a  magnetostrictive  actuator. 
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We  first  showed  that,  the  inversion  error  for  a  Preisach  operator-based  hysteresis  model,  either 
rate-independent  or  rate-dependent,  can  be  bounded  in  magnitude.  Furthermore,  the  bound  can  be 
characterized  in  terms  of  parameter  uncertainties  and  the  specific  inversion  scheme  used.  This  allows 
us  to  model  the  inversion  error  as  an  exogenous  noise  and  attenuate  its  impact  by  robust  control 
techniques. 

We  then  aimed  to  achieve  three  competing  goals  in  the  robust  control  problem  formulation:  ro¬ 
bust  stability,  robust  tracking,  and  not  exceeding  the  saturation  constraint.  We  also  discussed  how 
to  balance  these  goals  and  seek  the  sensible  trade-off.  Simulation  and  experimental  results  have 
demonstrated  the  effectiveness  of  the  robust  control  scheme. 
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